 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Loglinear model fitting in the poissonreg-proto
;;
;; Just say 
;; 
;; (log-linear-model ncat counts model vars)
;;
;; where ncat is a list with the number of categories of 
;; the variables, counts is a list of cell counts, model
;; is a list of lists indicating which interactions go
;; in the model, for example '((0 1) (2 3)) or '((0) (0 1)),
;; and the optional argument vars is a list of variable 
;; names, for example '(pre post effect).
;;
;; This file also provides the utility expand-hierarchy,
;; which expands a hierarchical model to a list of all
;; interactions. Thus
;; 
;; (expand-hierarchy '((0 1) (2 3)))
;;
;; returns
;;
;; ((0) (1) (2) (3) (0 1) (2 3))
;;
;; which can then be used as an argument to log-linear-model.
;;
;; In addition we use utilities for decoding and encoding
;; vectors, and for lexicographical ordering of lists of
;; indices. They are in the separate files decode-encode.lsp
;; and lexico-sort.lsp.
;;
;; Version 1.0, Barcelona, Spain, Jan de Leeuw, April 27, 1995
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;Additions Pedro Valero 2001

;; new methods for make-names and make-pred called my-make-names and
;; my-make-pred substituting original (apparently buggy) methods
;; Pedro Valero Feb 2002


(require "glim")
;(require "decode-encode")
;(require "lexico-sort")

(defproto log-linear-proto 
  '(model ncat b m) nil poissonreg-proto)

(defmeth log-linear-proto :isnew (the-ncat the-counts the-model the-vars 
                                           &optional level-labels 
                                           the-data-matrix
                                           reference-categories) ;PV added the argument data-matrix and reference-categories

  (send self :set-model the-model)
  (send self :set-ncat the-ncat)
(call-next-method 
 :y the-counts
 :x (send self :my-make-pred the-model the-data-matrix level-labels reference-categories)
 :predictor-names (send self :my-make-names the-model the-vars level-labels reference-categories)
))

(defmeth log-linear-proto :verbose (&optional verb)
  (if verb 
      t
      nil)
  )

(defmeth log-linear-proto :my-make-names (model vars level-labels &optional reference-categories)
  (let* 
    ((level-labels (if reference-categories 
                           (mapcar #'(lambda (cat label)
                                       (remove cat label :test #'equal)) 
                                   reference-categories level-labels)
                           (mapcar 'cdr level-labels)))
     (n-level (mapcar #'length level-labels))
     (step1 
      (mapcar #'(lambda (mod)  
                  (mapcar #'(lambda (var levels) 
                              (mapcar #'(lambda (variable)
                                          (mapcar #'(lambda (level) 
                                                      (strcat (princ-to-string variable) "(" (princ-to-string level) ")")) 
                                                  levels))
                                      (list var)))
                          (select vars mod) (select level-labels mod)))
              model))
     (step2 
      (combine 
             (mapcar #'(lambda (stp1 mod) 
                       (mapcar #'(lambda (prf)
                                   (introduce-symbol 
                                    (mapcar #'(lambda (ind pr) (select (combine ind) pr))
                                            stp1  prf)))
                               (all-profile-list (select n-level (combine mod)))))
 
                   step1 model)))
     )
    step2))
    
  


(defmeth log-linear-proto :my-make-pred (model data-matrix level-labels &optional reference-categories)
  "Compute matrix of indicators using the original data matrix in order to avoid errors for data no organized in a standard way. Removes first category"
  (let* ((level-labels (if reference-categories 
                           (mapcar #'(lambda (cat label) 
                                       (remove cat label :test #'equal))
                                   reference-categories level-labels)
                           (mapcar 'cdr level-labels)))
         (n-level (mapcar 'length level-labels))
         (nvars (length level-labels))
         (model model)
         ;(terms-in-model (length model))
         (data-matrix (mapcar #'(lambda (col)
                                  (coerce col 'list))
                                  (column-list data-matrix)))
         (n (length (first data-matrix)))
         (null-vector (repeat '0 n))
         (temp-vector)
         (indicators
          (mapcar #'(lambda (mod)
                      (mapcar #'(lambda (labels col)
                                  (mapcar #'(lambda (label) 
                                              (setf temp-vector (copy-list null-vector))
                                              (setf  (select  temp-vector
                                                              (which
                                                               (mapcar #'(lambda (val) (equal label val))
                                                                       col)))
                                                     1)
                                              temp-vector)
                                          labels)
                                  )
                              (select level-labels mod) (select data-matrix mod)
                              ))
                  model))
         )
    (apply 'append 
           (mapcar #'(lambda (indic mod) 
                       (mapcar #'(lambda (prf)
                                   (apply '* (mapcar #'(lambda (ind pr) (select ind pr))
                                                     indic   prf)))
                               (all-profile-list (select n-level (combine mod)))))
 
                   indicators model))
            
    ))
  


(defmeth log-linear-proto :make-pred (model)
  (let* ((nc (send self :set-ncat))
         (uu (all-profile-list nc))
         (xx nil))
    (dolist (mm model xx)
            (let* ((mn (select nc mm))
                   (vv (mapcar #'(lambda (x) (select x mm)) uu))
                   (v-values (remove-duplicates vv :test #'equal))
                   (tt (mapcar #'(lambda (vv-values)
                                   (if (member '0 vv-values)
                                       0
                                       (position vv-values v-values :test #'equal)))
                               vv)))
              ; (tt (mapcar #'(lambda (x) (decode x mn :zero-special t)) vv))) 31/8/2001 
            ;PV There is a bug with decode that produces wrong indicators.
           ; The line above works better
            (setf xx (append xx (indicators tt)))))
    ))

(defmeth log-linear-proto :make-names (model vars &optional level-labels) 
(let* ((nc (send self :set-ncat))
       (uu (all-profile-list nc))
       (xx nil)
       (level-labels level-labels))
(dolist (mm model xx)
(let* ((mn (select nc mm))
       (nn (iseq (length mm)))
       (vv (mapcar #'(lambda (x) (select x mm)) uu)))
  (when level-labels (setf this-level-labels (select level-labels mm)))
  (setf xx (append xx (expand-names mm vv vars this-level-labels)))
  

#|(setf xx (if level-labels 
             (append xx (expand-names mm vv vars level-labels))
             (append xx (expand-names mm vv vars))))|#
))))

(defmeth log-linear-proto :set-ncat (&optional (the-ncat nil set))
 (if set (setf (slot-value 'ncat) the-ncat)
     (slot-value 'ncat))
)

(defmeth log-linear-proto :set-model (&optional (the-model nil set))
 (if set (setf (slot-value 'model) the-model)
     (slot-value 'model))
)

(defmeth log-linear-proto :b (&optional b)
 (if b (setf (slot-value 'b) b )
     (slot-value 'b))
)

(defmeth log-linear-proto :m (&optional m)
 (if m (setf (slot-value 'm) m )
     (slot-value 'm))
)


;added by Pedro Valero 9/2000
(defmeth log-linear-proto :display ()
  )


(defmeth log-linear-proto :aic ()
  (let* ((dgres (- (apply '* (send self :set-ncat))
                 (send self :num-coefs)))
        (aic (- (send self :deviance) (* 2 dgres)))
        )
    aic)
  )

(defmeth log-linear-proto :bic ()
  (let* ((dgres (- (apply '* (send self :set-ncat))
                 (send self :num-coefs)))
        (bic (- (send self :deviance) (* dgres (log (sum (send self :yvar))))))
        )
    bic)
  )

(defmeth log-linear-proto :chi ()
  (sum (** (send self :chi-residuals) 2)))

(defmeth log-linear-proto :adj-deviance-residuals ()
  "Computes adjusted deviance residuals"
  (/ (send self :deviance-residuals) (sqrt (- 1 (send self :leverages))))
  )

(defmeth log-linear-proto :adj-chi-residuals ()
    "Computes adjusted Pearson residuals"
  (let (
         (leverages 
          (mapcar #'(lambda (val) (if (= val 1) .999999999 val))
                  (send self :leverages)))
         )

  (mapcar '(lambda (va) (if (complexp va) 0
                            va))
             (/ (send self :raw-residuals) 
     (sqrt (* (send self :raw-predicted)
              (- 1 leverages)))))
  ))


(defmeth log-linear-proto :change-chi ()
  (** (send self :adj-chi-residuals) 2))

(defmeth log-linear-proto :raw-predicted ()
  "Predicted values in original scale"
  (send self :fit-means)
  )

#|
This versions of leverages and cook distances give the same result as standard methods of glim and are not used. Just left to remind me the time wasted learning this.PV

(defmeth log-linear-proto :my-leverages ()
  "New leverages method for loglinear"
  (let ((x (send self :x-matrix))
        (w (diagonal (send self :weights)))
        (rtw (sqrt (diagonal (send self :weights)))))
                (diagonal (matmult rtw x
                                   (inverse (matmult 
                                             (transpose x) w x)) 
                                   (transpose x) rtw))))

(defmeth log-linear-proto :my-Cook-distances ()
  "Computed adjusted Cook distances for the loglinear model. "
  (let ((leverages (send self :leverages)))
 (/ (* (** (send self :adj-chi-residuals) 2) leverages)
    (* (send self :num-coefs) (- 1 leverages)))
    ))

|#
(defmeth log-linear-proto :leverages ()
  (let* (
         (temp (call-next-method))
         )
    (mapcar #'(lambda (x) (if (= x 1) 1.00000001 x))
            temp)))

(defmeth log-linear-proto :one-iter-uni-dim-newton ()
  (let* (
         (x (send self :x-matrix))
         (n (send self :yvar))
         (m-1 (if (send self :m) (coerce (send self :m) 'vector)
                  (repeat 1 (length n))))
         (b-1 (if (send self :b) (send self :b)    
                  (send (regression-model x
                                          (log (+ 0.5 (send self :yvar))) :display nil)
                        :coef-estimates)))
         #|(m-1 (if (send self :m) (send self :m) 
                  n))
         (b-1 (if (send self :b) (send self :b)    
                  (repeat 0 (length (column-list (send self :x-matrix))))))|#
                  
         (dif  (- n m-1))
         (den (coerce (matmult dif x) 'list))
         (coefs (mapcar #'(lambda (xi bj den)
                            (+  bj                                         
                                (log (+ 1 (/ den
                                             (sum (* m-1 (matmult xi xi))))))))
                        (column-list x) b-1 den))
         (m (exp (matmult x coefs))))
    (send self :b coefs)
    (send self :m m)
    (list coefs m)))

#|(defmeth log-linear-proto :correct-yvar ()
  (when (which (= (send self :yvar) 0))
       (send self :yvar
             (+ (send self :yvar) 0.5))))|#

(defmeth log-linear-proto :compute-uni-dim-newton ()
  (let* ((epsilon (send self :epsilon))
         (maxcount 200)
         (low-lim .0000001)
         (tempb)
         )
    (send self :one-iter-uni-dim-newton)
    (dotimes (i maxcount)
           (setf tempb  (send self :b))
         (send self :one-iter-uni-dim-newton)
             (when (< (sum (abs (- tempb (send self :b))))
                      low-lim)
                   (return))
           ;  (print i)
            ; (print (list i (max (abs (- tempb (send self :b))))))
             )))
  

(defun log-linear-model 
  (ncat counts model &optional (vars (iseq (length ncat)))  (level-labels nil) (data-matrix nil) (reference-categories nil))
  (send log-linear-proto :new ncat counts model vars level-labels data-matrix reference-categories))

(defun expand-hierarchy (model)
"Args: model
Takes a hierarchical model MODEL (list of lists) and
expands it to a full model (list of all unique
sublists of MODEL)"
  (let ((ex nil)
        (as (all-subsets (iseq (1+ (max model))))))
(dolist (mm as)
(if (is-subsetp mm model)
    (setf ex (append ex (list mm)))))
(lexico-sort (rest ex))
))

(defun expand-names (comp blow vars &optional label-levels)
(let ((n (length comp))
      (label-levels label-levels))
(if (= n 1)
    (setf blow (mapcar #'(lambda (x) (select (combine label-levels) x)) blow))
    (setf blow (mapcar #'(lambda (y)
                           (map-elements #'(lambda (x z)
                                             (select x z))
                                         label-levels y))
                       blow)))

(if (= n 1)
    (level-names (mapcar #'first blow) 
                 :prefix (elt vars (first comp)))
    (apply #'cross-names
           (mapcar #'(lambda (y) 
                       (level-names 
                        (mapcar #'(lambda (z) (elt z y)) blow)
                        :prefix (elt vars (elt comp y)))) (iseq n))))
))

;;;;;;;;;;;;;;;;;;follows decode-encode.lsp

;; Jan de Leeuw, April 27, 1995

(defun encode (x base)
"Args: number base
Convert decimal NUMBER to NUMBER in base BASE."
  (let* ((nb (length base))
         (lb (make-list nb)))
    (dotimes (i nb lb)
             (let* ((k (1- (- nb i)))
                    (b (elt base k))
                    (y (mod x b)))
               (setf (elt lb k) y)
               (setf x (/ (- x y) b))
               ))
 )) 

(defun decode (x base &key (zero-special nil))
"Args: number base
Converts base BASE number NUMBER to decimal number. Both
NUMBER and BASE are lists. If ZERO-SPECIAL is true, then
all numbers with a zero in it are transformed to zero."
(let* ((y (reverse (accumulate #'* base)))
       (z (rest (combine y 1))))
(if zero-special
    (if (member 0 x) 0 (sum (* x z)))
    (sum (* x z)))
));el fallo está aquí, hay que hacer que produzca números diferentes para cada celda. 

(defun all-profile-list (base)
"Args: base
Makes a list of lists, each list of the same length of base,
with all combinations of numbers between zero and base."
  (mapcar #'(lambda (x) (encode x base)) (iseq (prod base)))
)

(defun all-profile-matrix (base)
"Args: base
Makes an array, each row of the same length of base,
with all combinations of numbers between zero and base."
(let ((m (length base))
      (n (prod base)))
(make-array (list n m) :displaced-to
            (coerce (element-seq (all-profile-list base)) 'vector))
))

(defun all-subsets (list)
"Args: list
Power set of list."
(let* ((n (length list))
       (nn (^ 2 n)))
(mapcar #'(lambda (x) (select list (which (= 1 x))))
        (mapcar #'(lambda (x) (encode x (make-list n :initial-element 2)))
                (iseq nn)))
))

(defun is-subsetp (x y)
"Args: list1 list2
Is LIST1 a subset of one of the lists in LIST2 ?"
(nor (mapcar #'(lambda (z) (subsetp x z)) y))
)

(defun nor (list)
"Args: list
Returns true if at least one of the elements is true."
  (if (singlep list) (or (first list))
      (or (first list) (nor (rest list))))
)

(defun nand (list)
"Args: list
Returns true if all elements are true."
  (if (singlep list) (and (first list))
      (and (first list) (nand (rest list))))
)

(defun singlep (x)
"Args: x
Is X a list with a single element ?"
  (equal x (list (first x))))


;;;;;;;;;;follows lexico-sort

;; Jan de Leeuw, April 27, 1995

(defun lexico-sort (x)
"Args: x
Sorts a list of lists according to length, and
lists of the same length lexicographically with
respect to <."
  (sort x #'lexico-compare)
)

(defun lexico-compare (a b)
(cond 
  ((< (length a) (length b)) t)
  ((> (length a) (length b)) nil)
  (t (cond  
       ((< (first a) (first b)) t)
       ((> (first a) (first b)) nil)
       (t (lexico-compare (rest a) (rest b))))))
)